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Abstract 

The aim of this paper is to carry out a mathematical analysis of a 
system of ordinary differential equations introduced by R. Lev Bar-Or to 
model the interactions between T cells and macrophages. Under certain 
restrictions on the parameters of the model, theorems are proved about 
the number of stationary solutions and their stability. In some cases the 
existence of periodic solutions or heteroclinic cycles is ruled out. Evidence 
is presented that the same biological phenomena could be equally well 
described by a simpler model. 

1 Introduction 

Autoimmune diseases result in a great deal of suffering for affected individuals 
and huge costs for society. Notable examples are multiple sclerosis, rheumatoid 
arthritis and type I diabetes. In these diseases the ability of the immune system 
to distinguish between self and non-self is compromised, with the result that host 
tissues are attacked and damaged. It is important to get a better understanding 
of the processes involved and one way to do so is to introduce theoretical models 
of the immune system, in particular mathematical models. 

There has been a lot of work on mathematical modelling of interactions of the 
immune system with pathogens. See for instance the book of Nowak and May 
[7] which concentrates on the case of HIV. Autoimmune diseases do not need to 
involve any pathogens, although pathogens might contribute to them indirectly. 
For this reason it would be interesting to have models for the intrinsic workings 
of the immune system where non-self antigens play no direct role. Apparently 
few models of this type exist in the literature. One example is introduced in a 
paper of Lev Bar-Or [5j. It is a system of four ordinary differential equations 
which describes the interactions of T cells and macrophages by means of the 
cytokines they produce. The aim of this paper is to investigate what can be said 



about the properties of the solutions of this system on the level of mathematical 
proofs. 

T cells are white blood cells which mature in the thymus. One type of T 
cell, the T helper cell, helps to direct the activity of other immune cells. These 
cells are also known as CD4 + since they carry the surface molecule CD4. T 
cells communicate with other cells by secreting soluble substances known as cy- 
tokines. Another important type of blood cell is the macrophage which ingests 
pathogens and cell debris through phagocytosis. Macrophages also secrete cy- 
tokines. Both T cells and macrophages react in various ways to the cytokines 
which are present in their surroundings. The list of known cytokines is long 
and each of them has its own characteristics in terms of which types of cells 
secrete it and what effects it has on cells which detect its presence. An idea of 
the complexity of this signalling system can be obtained from [5]. 

It is common to distinguish between two types of T helper cells, known as 
Thl and Th2, according to the cytokines they produce. There may be overlaps 
but roughly speaking it may be supposed that there is one set of cytokines which 
are called type 1 and are produced by Thl cells and another called type 2 which 
are produced by Th2 cells. Macrophages produce cytokines of both types. The 
basic quantities in the equations of [5] are average concentrations corresponding 
to type 1 and type 2 cytokines produced by T cells and macrophages. These 
four concentrations are denoted by Cf , Cj, Cf 1 and C^ 1 ■ The populations of 
different cell types do not occur directly in the system. The quantities which 
have a chance of being measured directly are C\ — \{Cfy + Cf 1 ) and C% — 
\{C% + C?). 

There are cases where the immune response is dominated by either Thl or 
Th2 cells. This may be important in order to effectively combat a particular 
pathogen. For instance a sufficiently strong Thl response is necessary for con- 
taining or eliminating a tuberculosis infection [9] , [6] . An inappropriate balance 
between these two states can also contribute to autoimmune disorders. It has, 
for instance, been suggested that multiple sclerosis is associated with an immune 
response which is biased towards Thl. See [3] for a critical review of this idea. 
In the context of the model it is said that there is Thl or Th2 dominance if 
C\ > C2 or C2 > Ci, respectively. 

Another important function of macrophages which plays a role in the model 
of [5] is the presentation of antigens. Small peptides which result from the 
digestion of material taken up by a macrophage are presented on its surface in 
conjunction with MHC II molecules. (Major histocompatibility complex of class 
II.) This can stimulate T cells which come into contact with the macrophage. 

It has not been possible to give a complete analysis of the dynamics of 
the system of [5]. A certain inequality on the parameters of the system gives 
rise to a regime where there is a unique stationary solution which acts as an 
attractor for all solutions as t —> 00. This is proved in Theorem 1. For a more 
restricted set of parameters it is possible to show (Theorem 2) that each solution 
converges to a stationary solution which in general depends on the solution 
considered. With further restrictions on the parameters it is shown that there 
are between one and three stationary solutions and the subsets of parameters for 
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which different numbers of stationary solutions occur are described. This is the 
content of Theorem 3. In Theorems 2 and 3 the coefficients describing antigen 
presentation are set to zero. One situation where information can be obtained 
on the dynamics including antigen presentation is analysed in Theorem 4. 

The analysis which has been done has uncovered no evidence that the inclu- 
sion of the effect of macrophages makes an essential difference to the behaviour 
of solutions. The types of qualitative behaviour which have been proved to oc- 
cur in this paper and those which are shown in the figures in [5] can be found in 
a truncated system which only includes the effect of T cells. This is discussed 
in section [3] 

2 Analysis of the dynamical system 

In what follows it will be convenient to use a notation which is more concise 
than that of [5]. Let x\, X2, Xs, X4, z\ and denote Cj ', Cj, , C% , %i +£3 
and X2 + £4 respectively. The basic dynamical system is: 

dx ' 

-^ = -d i x i +g{h i ); i = l,2,3,4. (1) 
The di are positive constants. The function g is given by 

ff(s) = i(l + tanh(a:-0)) (2) 

where 9 is a constant. It satisfies the relations g(x + 6) + g{—x + 6) = 1 and 
g'(x) — 2g(x)(l — g(x)). Hence g{0) = h and g'{6) — h. The functions hi are 
defined by hi = Y]j aijXj for some constants a^-. The coefficients in (JTJ) satisfy 
the following conditions 

1. a,ij — bij + Cij, i — 1,2, for some coefficients bij and Cij 

2. b\j — —b2j for all j. 

3. 6n > 0, 613 > 0, 612 < and 614 < 

4. The ratio c%jjc\j is independent of j and positive. 

5. c\\ > and C13 > 

6. a^j — —ci4j for all j 

7. a 3 i > 0, a 33 > 0, a 32 < 0, a 34 < 

The sign conditions encode the fact that the effect of type 1 cytokines on cells is 
to increase their production of type 1 cytokines and to decrease their production 
of type 2 cytokines, while the effect of type 2 cytokines is exactly the opposite. 
The coefficients cy encode the effects of antigen presentation. No assumption 
is made on the signs of C12 and C14. There are eighteen parameters in the 
model which are only constrained by some positivity conditions. The quantities 
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z\ and Z2 represent total concentrations of type 1 and type 2 cytokines. The 
biologically relevant region B is that where all the Xi are non-negative. The 
function g is strictly positive. Hence if one of the variables Xi vanishes at some 
time its derivative at that time is strictly positive. It follows that B is positively 
invariant under the evolution defined by the system. If some Xi is greater than 
or equal to d^ 1 on some time interval then Xi is decreasing at a uniform rate 
during that time. It follows that all solutions exist globally to the future and 
enter the region B\ defined by the inequalities Xi < d^ 1 after finite time. Thus in 
order to study the late-time behaviour of any solution starting in B it is enough 
to consider solutions starting in B\. In fact any solution enters the interior of 
B\ after finite time. 

Theorem 1 For any value of the parameters the system (TTJ has at least one 
stationary solution. If 

sup \a%j\ < 2inf di. (3) 

i ■ % 
3 

then there is only one stationary solution and all solutions converge to it as 
t — » oo. 

Proof That the system always has at least one stationary solution follows from 
the Brouwer fixed point theorem, (cf. [I], Theorem 1.8.2). A stationary solution 
of the system satisfies 

Xi = d^gQii). (4) 
To prove the second part of the theorem consider the following estimate 

Id^gihiy^-d^gihiix))] < {2d i )- 1 \h l {y)-h l {x)\ < (2d j )~ 1 ^ 1^1% - Xj \. 

3 

(5) 

The first of these inequalities uses the mean value theorem and the fact the 
derivative of g is nowhere greater than one half. If ([3]) holds then the mapping 
{xi} i — ^ {d^ 1 g(hi)} maps B\ to itself and is a contraction in the maximum norm. 
Hence it has a unique fixed point. It follows that when the coefficients of the 
system satisfy the restriction ([3]) the system has exactly one stationary solution. 
If x(t) and y(t) are two solutions then under the assumption © it can be shown 
that | x — y\ decays exponentially as t — > oo. Thus all solutions converge to the 
unique stationary solution as t — > oo. A related statement is that if cc* is the 
unique stationary solution then | Lyapunov function. 

A limiting case of ([T]) is that where all vanish. This will be called the zero 
MHC system. The pattern of signs in the coefficients in the zero MHC system is 
such that the change of variables Xi — (— l) t+1 Xi leads to a system ^ = f(xj) 
satisfying > for all i ^ j. This means that the dynamical system for 
the Xi is cooperative [2] . This pattern of signs corresponds to what is referred 
to as a 'community with limited competition' or 'competing subcommunities 
of mutualists' in [5]. It implies, using the Perron- Frobenius theorem, that the 
linearization at any point of the vector field defining (fTJ) has a real eigenvalue 
of multiplicity one which is greater than the modulus of any other eigenvalue. 
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Notice that in the special case of the zero MHC system where all di are equal 
to d the linearization has two eigenvalues equal to —d. As a consequence all 
the eigenvalues of the linearization must be real in this case. It is not clear how 
these facts about the eigenvalue structure can be used to help to understand 
the dynamics. The condition on the coefficients di, which says that the rate 
of degradation of different cytokines is exactly equal, would not be true with 
biologically motivated parameters. Nevertheless it is not unreasonable to assume 
that these coefficients are approximately equal and that the model with equal 
coefficients is not a bad approximation to the situation to be modelled. An 
assumption of this type is made in the model of [5]. 

Call the system obtained by assuming = 0, d\ = d%, d 3 = d^ and 9 = 
in ((T|) System 2 while (JXJ) itself is System 1. The condition 8 = implies that 
g{x) + g{—x) = 1. Adding the equations for i = 1 and i = 2 then gives 

— (xi + x 2 ) = -di(xi +ar 2 ) + 1. (6) 

Similarly 

— (x 3 +X4) = -d 3 (x 3 +Xi) + 1. (7) 

These equations can easily be analysed. An invariant manifold Si is defined by 
x\ + X2 = d^ 1 and x 3 + X4 — d^ 1 . Substituting this back into the full system 
gives a two-dimensional system written out below which will be called System 3. 
The w-limit set of any solution of System 2 is contained in the invariant manifold 
Si . Passing from System I to System 3 means setting some parameters to zero 
and then restricting to a two-dimensional invariant manifold of the resulting 
system. System 3 may not include all the most interesting dynamics exhibited 
by solutions of the system of [5] but can be perturbed to get information about 
cases where the effect of the MHC is small but non-zero. The phase portrait 
of System 2 follows immediately from that of System 3. The explicit form of 
System 3 is: 

dixi + g((au ~ ai 2 )xi + (a i3 - a u )x 3 + a^d^ 1 + aud^ 1 ) (8) 
^3^3 + 5(( a 3i ~ a-32)xi + (a 33 - 034)2:3 + a^d^ 1 + a^d^ 1 ) (9) 

Note that the coefficients of xi and X3 in these equations which are linear 
combinations of the ay are all positive. The constant terms are negative. This 
is a cooperative system. This fact together with the fact that the dimension of 
the system is two implies that each solution converges to a limit as t —> 00 [2] . In 
other words the w-limit set of each solution is a single point. Furthermore, there 
are no homoclinic orbits or heteroclinic cycles. What is not clear in general is 
how many stationary solutions there are. Some of the conclusions of the above 
discussion can be summarized as follows: 

Theorem 2 If cy = for all i, j , di = d 3 — d^ and = then every 
solution of Q]) converges to a stationary solution satisfying x\ + x 2 — d^ 1 and 
x 3 + X4 = d^ 1 as t — > 00. There are no homoclinic orbits or heteroclinic cycles. 



dxi 

~dt 

dx 3 

~~dt 
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Consider now the special case of System 3 obtained by setting ay = 0,3 j 
for all j and d\ = d 3 . Call it System 4. This corresponds to the assumptions 
that the T cells and macrophages have identical properties with respect to their 
death rate and their interactions with cytokines. This system consists of two 
copies of a single equation. The quantity x\ — x% decays exponentially. The 
further assumption x\ = £3 on the initial data, which means that there are 
equal numbers of T cells and macrophages, gives rise to a single ODE, call it 
the toy model. It is the same equation which occurs twice in System 4. It is of 
the form 

dx 

— = —d%x + g(ax — 6) (10) 

where a = an — ai2 + 013 — an and 6 = — (ai2 + ai^dj -1 . Note that a and b are 
positive. Here the zero MHC condition has been used. This can be simplified 
by defining x' = d\X, t' = d\t and a' = Suppressing the primes leads to 

-77 = -x + g(ax - b) (11) 

Denote the right hand side of (jlll) by h(x), suppressing the parameter depen- 
dence. Stationary solutions of the toy model are given by solutions of the 
equation 

g(ax — b) = x. (12) 

Since < g < 1 all solutions of this equation are contained in the interval (0, 1). 
The right hand side of this equation takes the value zero at x = and the value 
one at x = 1. Thus by the intermediate value theorem (fT2j) has at least one 
solution on the interval of interest. Now 

— (g(ax — b) — x) = 2ag(ax — 6)[1 — g(ax — 6)1 — 1. (13) 
dx 

The first term on the right hand side of (|13[) is symmetric about x = — where it 
attains its maximum value § . Thus this derivative vanishes at zero, one or two 
points of the real line according to whether a < 2, a = 2 or a > 2. For a > 2 call 
the zeroes of the derivative x\ and X2 with x\ < xi. It can be concluded that 
for any value of a there are at most three solutions of (|12l) . All these solutions 
must be within the interval of interest. If there are exactly two solutions then 
one of them must correspond to a point where the right hand side of (|11[) and its 
derivative vanish simultaneously. For otherwise the value of this function at one 
would be negative for large x, a contradiction. For a given value of a > 2 there is 
precisely one value of b for which the equation —x\ -I- g(ax% — b) = holds. Call 
it 61(a). Define 6 2 (a) similarly, replacing x% by Then bi(a) and 62(a) are the 
only values of 6 for which there is a simultaneous solution of g(ax — 6) = x and 
ag'{ax — 6) = 1. For all a > 2 the inequality 61(a) < 62(a) holds. As a —> 2 both 
tend to one. The union of the graphs of 61 and 62 divides the (a, 6)-plane into 
two regions. On each of these the number of stationary solutions is constant 
and equal to one or three. On the region including points with a < 2 it is 
obviously one. Considering points where 6 = \ and a slightly larger than two 
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makes makes it clear that it is three on the other region. Information can also 
be obtained by simple geometric considerations on the position of the stationary 
points. If a < 2 and 6 < ~ then the stationary point has x < i while for 6 > | 
it satisfies x > |. If a > 2 and there are three stationary points then the central 
(unstable) one satisfies x < ^ for b < | and x > | for b > §. 

When there is only one solution it is a hyperbolic sink. When there are three 
solutions the two outermost are hyperbolic sinks while the intermediate one is 
a hyperbolic source. If the parameter b is held fixed at some value less than one 
and a is varied then at the value of a where b\ (a) = b there is a fold bifurcation. 
To see this note that 

-l + ag'(ax-b) (14) 

xg' (ax — b) = 2xg(ax — 6)(1 — g(ax — b)) (15) 

= a 2 g"(ax -b)= 4a 2 g'(ax - 6)(1 - 2g(ax - b)) (16) 

From the first equation it follows that at a critical point x the quanitity ag'(ax — 
b) must be equal to one. Hence the derivative with respect to a does not 
vanish there. The second derivative with respect to x can only vanish there if 
g(ax — b) = |, which implies that x — -, x = ^ and 6=1. The situation for 
6 > 1 is similar to that for 6 < 1 . To understand the case 6=1 note that there 
is a cusp bifurcation in (a, 6) at the point (2, 1). This can be shown using the 
facts that 

-g'(ax-b), (17) 
~ag"{ax - b), (18) 
axg" (ax — 6), +g'(ax — 6) 

g"'(ax - 6) (19) 
4a 3 g"(ax - 6)(1 - 2g(ax - 6)) - 8a 3 (g'(ax - 6)) 2 . (20) 

Hence §^4 and ^-§~ik — ^k^-s- are both non- vanishing and Theorem 8.1 of G2 

ar oa axab ob oxoa L - 

applies. 

If we remove the condition X\ — xj, then x\ — 2:3 decays exponentially. The 
qualitative nature of the dynamics of System 4 with the given restrictions on 
the parameters is then clear. This in turn gives information about the phase 
portrait of System 2 with these values of the parameters. Converting back to 
the original variables has the effect of replacing a by ajdi in these criteria. Some 
of these results are summarized in the following theorem: 

Theorem 3 If in Theorem 2 it is additionally assumed that aij = a%j for all j 
and d\ = (I3 then the number of stationary solutions of ([T]) is between one and 
three for any values of the parameters. It is one whenever an— <Zi2+<2i3— &14 _ 2. 



dh 

dx 
dh 

da 

d 2 h 

dx 2 



dh 
db 
d 2 h 
dxdb 
d 2 h 
dxda 
d 3 h 
dx 3 
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There is a non-empty open set where it is three. The set where it is two is a 
union of two smoothly embedded curves which is the boundary between the sets 
where it is one and three. 

Consider a choice of parameter set for ([T]) for which there are three hyperbolic 
stationary points. It follows by a simple stability argument that there is an 
open neighbourhood of this parameter set in parameter space such that for all 
parameters in this neighbourhood there are precisely three hyperbolic stationary 
points. If the original set of stationary points consists of two sinks and one saddle 
then there is an open neighbourhood where this property persists. 

In the context of this theorem the question of Thl or Th2 dominance can 
be examined for stationary solutions. It is of interest to know whether changing 
the parameters in the system can cause a switch from one type of dominance 
to the other. If this happens there must be some values of the parameters for 
which there is a stationary solution with z\ = Zi- Under the hypotheses of the 
Theorem 1x2 — d^ 1 — x\ and x 4 = d^f 1 — x 3 for any stationary solution. Hence 
Z2 = 2e?^ 1 — z\. If z\ = Z2 then this implies that z\ = d^ 1 and x% — \d^[ . In 
terms of the new variable introduced in (fTTj) this means that x = 4. The above 
discussion gives information about Thl or Th2 dominance for various parameter 
values satifying the restrictions in the statement of Theorem 3. What is of most 
interest for the applications is the position of the stable stationary points. In 
particular it is clear that there exist parameter values where there are stationary 
solutions with both types of dominance. 

Now another set of simplified versions of the system will be studied. The 
phase portraits given in Fig. 2 in [5] relate to one special case of this kind. In the 
notation used here it is defined by the following restrictions on the parameters: 
di = 1 for all i, 8 = 0, all coefficients fry with j = 1 or j = 3 are equal, all 
coefficients a,-j with j = 2 or j = 4 are equal, the coefficients Cy are equal for 
all i and j. This reduces the number of free parameters to three. To simplify 
the notation let A — bn, B = —612, C = c\\. These are all positive constants. 
For Fig. 2 in [5] two sets of values for the coefficients are considered. In both 
cases C = 0.5. In Fig. 2a A = 0.4 and B = 0.5 while in Fig. 2b A = 0.6 and 
B = 0.65. With these assumptions System 1 becomes: 



dx\ 
~~dt 
dx 2 
~~dt 

dx 3 

dt 
dxi 

~dt 



-xi + g{{A + C)(x 1 + x 3 ) + (-B + C){x 2 + x A )) (21) 

-x 2 + g((-A + C)(xi + x 3 ) + (B + C)(x 2 + a*)) (22) 

-x 3 + g(A(x 1 +x 3 )-B{x 2 +x 4 )) (23) 

-x 4 + g(-A(x 1 +x 3 )+B{x 2 +x 4 )) (24) 



Call this System 5. If A + B + 2C < 1 then there is at most one stationary point 
of this system. This can be proved in the same way as Theorem 1. Adding the 
equations in pairs leads to a closed system for the two experimentally accessible 
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quantities Z\ and z 2 . These are the variables which are plotted in Fig. 2 of [5]. 



~dt 

dz 2 

dt 



zi + g((A + C)zx + (-B + C)z 2 ) + g(Azx - Bz 2 ) 



z 2 + g((~A + C)z x + {B + C)z 2 ) + g(-A Zl + Bz 2 ) (26) 



(25) 



Call this System 6. Denote the right hand sides of (|25|) and (|26l) by fx {z\,z%) 
and f 2 {z\, z 2 ) respectively. For this system uniqueness of the stationary solution 
follows from the assumption that the coefficients satisfy A + B + C < 1 . If in 
addition the coefficient C is assumed to vanish then this reduces to 



Call this System 7. Note that the system obtained by setting C = in System 
5 is equivalent to a special case of System 2 and that Theorems 2 and 3 apply 
to it. The parameters are related by a = 2{A + B), b = 2B. The special cases 
in Fig. 2 of [5] correspond to a = 1.8, b = 1 and a = 2.5, b = 1.3. It follows 
from ([27]) and (J2HJ) that f t (z x + z 2 ) = -{z\ + z 2 ) + 2 and that z\ + z 2 tends 
to two as t — > oo. Thus the dynamics is controlled by that on the invariant 
manifold Z\ = 2 — z 2 . With the choices which have been made it coincides with 
the invariant manifold S\ introduced earlier. 

Now consider what happens when C^O. Provided C < min{^4, B} the sys- 
tem consisting of (|25|) and (f26| is competitive. Thus each solution (zi, z 2 ) must 
converge to a stationary solution as t — > oo [2] . It can be concluded that the cor- 
responding solution {xi,X2,x$,Xi) tends to a stationary solution. Information 
about stationary points of the system can be obtained by stability considera- 
tions, as was mentioned following Theorem 3. This discussion is summed up in 
the following theorem: 

Theorem 4 Any solution Xi (i = 1,2,3,4) of ([2"Tj) -(124 |) with parameter values 
satisfying C < min{A, B} converges to a stationary solution as t — > oo. For 
fixed values of A and B and C sufficiently small the system has at most three 
stationary solutions and at most two stable stationary solutions. 

The results which have been obtained are unfortunately not sufficient to give a 
rigorous confirmation of the qualitative behaviour shown in Fig. 2a and Fig. 2b 
of [5]. The condition C < min{A,B} is not satisfied by the parameter values 
in Fig. 2a. It is satisfied in the case of Fig. 2b of [5] but no information is 
obtained about the number of solutions. 

3 Dynamics in the absence of macrophages 

In this section it will be shown that all the dynamical behaviour which has 
been shown to occur in system |T]) also occurs in a truncated model where the 



dz\ 
~dt 
dz 2 
~dT 



zi + 2g{Az x - Bz 2 ) 



z 2 + 2g{-Az x + Bz 2 ) 



(27) 



(28) 
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influence of macrophages is ignored. The system has two unknowns x\ and x 2 
with the same interpretation as before. The effect of the macrophages is turned 
off by setting et.y = when i is one or two and j is three or four. This leads to 
the following closed system for X\ and x 2 : 

d\x\ + g(a n xi - a 22 x 2 ) (29) 

d 2 x 2 + g(-anxi + a 22 x 2 ) (30) 

All the coefRents di, d 2l an and a 22 are assumed positive, as before. Restricting 
further to the case that di = d 2 = d and 6 = allows the analysis of the 
dynamics to be reduced to that of a scalar equation as in the previous section. 
The scalar equation is 

= -dxi + g((au + a 22 )xi - d~ 1 a 22 ). (31) 

This is essentially the equation (fTTj) analysed before and so all the phenomena 
found previously occur here also. 



ax i 

~dt 

dx 2 

~~dt 



4 Further remarks 

In this paper it has been possible to give a rigorous analysis of the asymptotics 
of solutions of the system of Lev Bar-Or [S] for some values of the parameters. 
Unfortunately there are large ranges of the parameters for which no conclusions 
were obtained or for which those conclusions are incomplete. The latter state- 
ment even includes the two cases for which numerical plots were included in [5] . 
The following questions have not been answered for the system ((T|) with general 
parameters. 

• Are there periodic solutions? 

• Are there homoclinic orbits? 

• Are there heteroclinic cycles? 

• Are there strange attractors? 

It should be emphasized that there is no evidence, analytical or numerical, that 
the answer to any of these questions is yes in the general case. We are left 
with a picture of a dynamical system where the long-time behaviour seems to 
be simple but it is quite unclear how to prove it except for a restricted set of 
values of the parameters. 

It is consistent with everything which has been found here that the system of 
four equations with eighteen parameters produces no qualitatively new phenom- 
ena in comparison with a reduced two-dimensional system with three param- 
eters. No effects were found which are specifically dependent on the inclusion 
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of the presentation of antigen by macrophages. The interactions between Thl 
and Th2 cells appear to be sufficient. It is seen that under some circumstances 
the model with T cells and cytokines alone predicts a situation of bistability 
where either a Thl or Th2 dominated state can be approached, depending on 
the initial data. 
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